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SYMPLECTIC  ATTITUDE  ESTIMATION  FOR 
SMALL  SATELLITES 


James  M.  ValpianD  &  Philip  L.  Palmer* 


In  this  paper,  a  novel  method  for  efficient  high-accuracy  satellite  attitude 
estimation  is  presented.  Symplectic  numerical  methods  are  applied  to  the  Ex¬ 
tended  Kalman  Filter  (EKF)  algorithm  to  give  the  SKF,  which  outperforms 
the  standard  EKF  in  the  presence  of  nonlinearity  and  low  measurement  noise 
in  the  1-D  case.  Building  on  this  result,  a  six-state  SKF  is  compared  to  an 
EKF  of  the  same  order  for  satellite  attitude  estimation.  Simulation  of  a  stan¬ 
dard  small  satellite  mission  demonstrates  orders  of  magnitude  improvement 
in  state  accuracy  and  preservation  of  constants  of  motion.  This  new  method 
shows  promise  for  improved  attitude  estimation  onboard  resource-constrained 
small  satellites. 

INTRODUCTION 

Increasing  satellite  attitude  requirements  demand  high  accuracy  estimation  methods 
capable  of  operating  under  significant  computational  constraints.  To  meet  these 
demands,  dynamical  modeling  has  been  used  as  an  effective  alternative  to  costly 
and  resource  intensive  ADCS  hardware  for  small  satellite  missions.  However,  these 
models  have  generally  been  propagated  using  methods  ill-suited  to  take  advantage 
of  the  unique  properties  of  Hamiltonian  systems. 

Recently  developed  symplectic  numerical  methods  have  demonstrated  promising 
results  when  applied  to  satellite  attitude  propagation.  Symplectic  methods  inher¬ 
ently  preserve  the  energy  of  Hamiltonian  systems  and  as  a  result  they  have  demon¬ 
strated  increased  stability  and  improved  long-term  performance  in  comparison  to 
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nonsymplectic  methods  (such  as  Runge-Kutta).  The  result  is  a  dramatic  improve¬ 
ment  in  state  propagation  accuracy  and  constants  of  motion  preservation  even  for 
symplectic  methods  inferior  in  order  to  nonsymplectic  ones.  Consequently,  they  are 
well  suited  for  systems  where  preservation  of  conserved  quantities  is  important  for 
accurate  predictions. 

Motivated  by  a  need  for  highly  accurate  and  efficient  attitude  estimation  al¬ 
gorithms,  this  paper  presents  a  novel  solution  that  capitalizes  on  the  strengths  of 
symplectic  methods.  Beginning  with  the  1-D  case,  a  symplectic  propagator  is  com¬ 
bined  with  Extended  Kalman  Filter  (EKF)  equations  to  give  a  symplectic  Kalman 
Filter  (SKF).  The  SKF’s  improved  performance  over  the  EKF  in  the  presence  of 
nonlinearity  and  low  measurement  noise  motivates  the  development  of  a  six-state 
SKF  algorithm  for  satellite  attitude  estimation.  Comparisons  with  a  standard  EKF 
demonstrate  significantly  better  steady-state  performance. 

SYMPLECTIC  NUMERICAL  METHODS 

A  symplectic  integrator  is  the  exact  solution  to  a  discrete  perturbed  Hamiltonian 
system  that  is  close  to  the  actual  Hamiltonian  of  interest.  Because  Hamiltonian 
systems  have  volume-preserving  mappings  in  phase  space  as  a  consequence  of  Liou- 
ville’s  Theorem,  symplectic  integrators  also  preserve  volume  in  phase  space  to  within 
machine  precision.  Therefore,  trajectories  modeled  by  these  integrators  do  not  cross 
in  phase  space  and  energy  errors  about  the  true  energy  are  bounded.  Practically 
speaking,  symplectic  integrators  preserve  the  invariants  of  motion  for  Hamiltonian 
systems  to  a  much  higher  degree  of  accuracy  than  non-symplectic  methods  of  the 
same  order.  When  applied  to  satellite  motion  which  is  governed  by  rigid  body 
dynamics,  symplectic  methods  exactly  conserve  angular  momentum  and  rotational 
energy  with  low  amplitude  sinusoidal  variation  and  no  numerical  dissipation. 

Generally,  non-symplectic  methods  will  introduce  secular  errors  into  conserved 
values  with  each  evaluation  and  therefore  misrepresent  the  conservative  system  as 
a  dissipative  one.  Over  time,  the  effect  of  this  dissipation  can  significantly  degrade 
state  accuracy.  In  contrast,  symplectic  integrators’  unique  properties  enable  them 
to  achieve  state  accuracies  equivalent  to  nonsymplectic  integrators  while  requiring 
less  function  evaluations  (i.e.  larger  time  steps)  to  do  so.  This  makes  them  ideally 
suited  to  model  Hamiltonian  systems  such  as  satellite  attitude,  particularly  onboard 
satellites  with  constrained  computational  resources.  When  computational  resources 
are  not  a  significant  concern,  the  comparatively  higher-accuracy  solutions  of  sym¬ 
plectic  integrators  reduce  the  requirement  for  frequent  measurement  updates  in  the 
attitude  estimation  problem. 

Recently,  symplectic  methods  have  been  applied  in  a  wide  variety  of  fields,  rang¬ 
ing  from  chemistry  to  celestial  mechanics.1  However,  these  methods  have  generally 
been  limited  to  stochastic  applications  where  the  behavior  of  a  multi-body  system 
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is  of  greater  concern  than  the  behavior  of  any  particular  system  component.2  Con¬ 
ventional  wisdom  held  that  symplectic  methods  would  not  be  useful  for  modeling 
specific  trajectories.3,4 

This  was  dispelled  by  Palmer  et  al5  in  1999,  when  they  published  an  orbit  estima¬ 
tion  algorithm  that  used  symplectic  methods  for  modeling  individual  satellite  orbits. 
Motivated  by  the  relatively  high  power  consumption  of  GPS  receivers  onboard  small 
satellites,  the  authors  developed  a  method  that  was  both  fast  and  highly  accurate 
to  minimize  CPU  power  consumption  and  reduce  the  number  of  GPS  measurements 
required  for  orbit  determination.  Using  symplectic  methods,  they  demonstrated 
modeling  accuracies  on  the  order  of  centimeters  for  propagations  spanning  multiple 
days  with  sparse  measurements. 

The  application  of  symplectic  numerical  methods  to  satellite  ADCS  design  is  a 
very  recent  development.  In  2004,  the  same  authors  presented  results  from  the  first 
application  of  symplectic  methods  to  satellite  attitude  propagation.6  The  authors 
sought  to  address  the  drawbacks  incurred  from  linearization  inherent  in  the  EKF 
approach  to  small  satellite  attitude  estimation,  including:  shorter  timesteps  and 
increased  computational  demand  needed  to  limit  error  growth  in  simplified  state 
prediction  methods;  limited  ability  to  utilize  more  accurate  integration  schemes  due 
to  constrained  computer  resources;  and  increased  polling  of  attitude  sensors  in  order 
to  compensate.  Noting  the  predictability  of  small  satellite  dynamics  and  the  small 
disturbances  they  experience,  the  authors  pointed  out  that  symplectic  methods  are 
well-suited  for  attitude  propagation.  Using  a  simple  low-order  composite  symplectic 
method  and  a  time-transformation,  they  propagated  torque-free  satellite  motion 
while  conserving  integrals  of  motion  exactly.  In  the  presence  of  gravity-gradient 
torque,  the  authors  demonstrated  low  amplitude  sinusoidal  errors  about  conserved 
system  values  using  large  timesteps;  in  contrast,  a  standard  nonsymplectic  Runge- 
Kutta  method  gave  secularly  increasing  errors  of  the  same  order  of  magnitude  using 
a  stepsize  two  orders  of  magnitude  smaller. 


1-D  SYMPLECTIC  KALMAN  FILTER 

Based  on  the  advantages  of  using  symplectic  numerical  methods  outlined  above,  it 
seemed  reasonable  to  hypothesize  that  they  provide  improved  accuracy  when  applied 
to  nonlinear  estimation  for  Hamiltonian  systems.  In  order  to  test  this  hypothesis 
and  determine  the  advantages  and  disadvantages  of  this  approach,  the  application 
of  symplectic  propagation  to  nonlinear  estimation  in  the  1-D  case  was  investigated. 
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1-D  EQUATIONS  OF  MOTION 

Starting  with  the  equations  for  a  nonlinear  pendulum 


H  —  ^p2  —  ft2  cos  (q) 

(la) 

Q  =  P 

(lb) 

p  =  — fi2  sin  ( q ) 

(lc) 

a  symplectic  discretization  was  applied.  The  numerical  method  chosen  was  the 
second-order  implicit  midpoint  rule  (IMPR)2 

Zj+i  =  Zj  +  Atf  (2) 

where  f  (z)  =  4|.  Note  that  subscripts  j  are  used  to  indicate  propagation  times. 
Applying  the  IMPR  method  to  Equations  (1)  gives 

(3a) 
(3b) 

(3c) 
(3d) 

In  this  formulation,  all  system  nonlinearity  is  contained  in  the  variable  Xi+i  and 
the  update  equations  for  qJ+i  and  pJ+i  are  simple  linear  equations.  The  iterative 
Newton-Raphson  method  was  employed  to  solve  Equation  (3d)  to  within  a  user- 
defined  tolerance. 


Qj+ 1  —  2Xj- (_i  qj 

Pj+i  =  pj  —  A tfl2  sin  (Aj+i) 
_  Qj  +  Qj+i 
2 

At 


Xj+i  = 


■ 

Xj+i  +  — H2  sin  (JCj+i)  =  ^Pj  +  Qj 


1-D  SYMPLECTIC  KALMAN  FILTER  EQUATIONS 

The  well-known  Extended  Kalman  Filter  algorithm  is  summarized  below.  Note  that 
the  subscripts  n  are  used  to  indicate  observation  times: 


•  Propagate  the  previous  time’s  state  estimate  x„_i  to  obtain  the  current  time’s 
state  prediction  x~ 


x~=x„_i  +  /  /  (xn_i)  dt 

Jtn~  i 


(4) 


•  Propagate  the  previous  time’s  error  covariance  estimate  Pn-i  to  obtain  the 
current  time’s  error  covariance  prediction  P~  using  the  state  transition  matrix 
<f>n  and  process  noise  Qn 

P-  =  $nPn~l$l  +  Qn  (5) 
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•  Calculate  the  Kalman  gain  Kn  using  the  observation  matrix  Hn  and  the  ob¬ 
servation  noise  covariance  Rn 

Kn  =  P~Hl  {Hnp-Hl  +  Rny 1  (6) 

•  Calculate  the  current  state  estimate  x„ 

x„  =  x“  +  Kn  (z„  -  Hnk~)  (7) 

•  Calculate  the  current  error  covariance  estimate  Pn 

Pn  =  (l-KnHn)p-  (8) 

The  standard  EKF  was  implemented  using  a  second-order  nonsymplectic  Runge- 
Kutta  numerical  method7  for  state  propagation  in  Equation  (4).  The  symplectic 
EKF  (hereafter  referred  to  as  SKF)  was  created  using  Equations  (3)  for  state  prop¬ 
agation. 

MONTE  CARLO  SIMULATION  SETUP 

In  order  to  compare  the  performance  of  the  EKF  and  the  SKF,  a  Monte  Carlo 
simulation  was  developed  using  the  following  conditions: 

•  State  vector  composed  of  generalized  coordinate,  q,  and  generalized  momen¬ 
tum,  p 

x  =  [q  pf  .  (9) 

•  Direct  observation  of  q 

H=[  10]  (10) 

•  Approximate  integration  of  the  state  transition  matrix 

dq  dq 

*  =  I2*2+  |  |  ^  (11) 

dq  dp 

•  Following  the  standard  practice  for  comparing  estimators  set  by  Athans,8,9,10 
process  noise  was  omitted  to  prevent  masking  of  linearization  errors 

Q  —  02x2  (12) 

To  determine  the  relative  merits  of  the  SKF  versus  the  EKF,  the  simulation  varied 

over  a  set  of  parameters  as  outlined  in  Table  1. 
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Nonlinearity  (q0,  rad) 

Measurement  Noise  Variance  (R,  rad2) 

Measurements/Orbit 

1.29 

(5  *  10-5)2 

1.86 

(5  *  10-4)2 

2.44 

(5  *  10~3)2 

100 

2.73 

(5  *  10~2)2 

50 

3.14 

(5  *  10-1)2 

10 

Table  1:  1-D  Monte  Carlo  Simulation  Parameter  Space 


The  expected  performance  of  both  filters  decreased  down  each  column  in  Table 

1,  with  the  pathological  parameter  set  defined  by  the  bottom  row.  For  each  of 
the  125  possible  combinations  of  parameters,  25  Monte  Carlo  runs  were  conducted. 
Each  run  used  unique  observation  sets  defined  by  Zrun  =  {zn  :  t0  <  tn  <  tend} 
where  zn  =  qn  +  vn  and  vn  ~  N  (0,  Rparameter  set)  ■  In  addition,  parameters  that  were 
constant  for  both  the  EKF  and  SKF  throughout  the  simulation  are  outlined  in  Table 

2. 


Propagation  At  (sec) 

Truth  At  (sec) 

jgjgyg 

i 

D2 

IMPR  Tolerance  (rad) 

5  *  10-3 

5  *  10“5 

1 

EH 

■ 

D 

1  *  HT12 

Table  2:  1-D  Monte  Carlo  Simulation  Constant  Parameters 


Truth  At  was  the  propagation  step  used  for  determining  the  true  values  of  the 
pendulum  state  and  IMPR  Tolerance  was  the  value  used  in  the  Newton- Raphson 
method  for  Equation  (3d).  Since  only  steady-state  performance  was  evaluated,  P0 
was  kept  constant  and  artificially  large  to  ensure  initial  filter  convergence. 

MONTE  CARLO  SIMULATION  RESULTS 

For  each  parameter  set,  the  absolute  value  of  the  position  error  was  averaged  across 
the  25  runs  at  each  observation  time.  The  results  from  the  filters  within  each  pa¬ 
rameter  set  fell  into  one  of  four  categories:  the  filter  converged,  the  filter  diverged, 
the  filter  failed  to  reach  steady-state  and  had  unbounded  error  growth  at  the  end  of 
the  simulation,  or  the  filter  reached  steady-state  with  errors  greater  than  3  a  of  the 
measurement  noise.  For  this  system,  99.73%  of  the  observations  are  within  3  a  of  the 
true  value;11  subsequently,  results  with  steady-state  errors  greater  than  this  value, 
results  with  unbounded  error  growth,  and  divergent  results  were  grouped  together 
as  undesirable  filter  results.  Then,  a  visualization  tool  was  used  to  make  general 
trend  observations  regarding  the  performance  of  each  filter  as  different  parameters 
varied;  results  are  shown  in  Figure  1.  These  visualizations  are  useful  for  developing 
an  intuitive  feel  for  the  filters’  relative  performance  trends.  Comparing  the  two  vi¬ 
sualizations,  it  appears  that  the  SKF  gives  convergent  behavior  in  regions  of  high 
nonlinearity,  high  measurement  frequency,  and  low  measurement  noise  whereas  the 
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1-D  SKF  Performance  Trends 


1-D  EKF  Performance  Trends 


Figure  1:  Comparison  of  1-D  SKF  and  EKF  Performance  Trends 


EKF  appears  to  fall  short  in  these  areas.  In  the  EKF  region  of  undesirable  perfor¬ 
mance  (red  volume) ,  the  ratio  of  the  position  error  below  which  the  filters  performed 
for  95%  of  the  steady-state  gave  a  relative  measure  of  the  advantage  of  using  the 
SKF  in  lieu  of  the  EKF.  For  example,  Figure  2  shows  the  effect  of  reducing  the 
measurement  noise  standard  deviation  by  an  order  of  magnitude  from  the  left  plot 
to  the  right:  the  ratio  of  the  95%  steady-state  error  value  increases  nearly  tenfold. 
Though  these  ratios  are  meaningless  as  independent  measures  of  performance  (due 
to  the  divergent  nature  of  the  EKF  solution),  they  demonstrate  the  advantage  of 
using  the  SKF  as  measurement  noise  decreases.  Similar  relationships  hold  for  in¬ 
creasing  nonlinearity  and  increasing  measurement  frequency;  these  results  indicate 
that  the  SKF  is  more  suited  than  the  EKF  to  handle  nonlinear  Hamiltonian  systems 
in  the  presence  of  high-accuracy,  high-frequency  observations. 


Figure  2:  Effect  of  Parameter  Variation  on  Relative  SKF  Advantage 


7 


SYMPLECTIC  KALMAN  FILTER  FOR  SATELLITE  ATTITUDE 


Building  on  the  results  from  the  1-D  case,  a  symplectic  propagator  was  developed 
based  on  previous  work  by  Palmer6, 12, 13  and  incorporated  into  a  six-state  SKF. 
Then,  comparisons  to  a  nonsymplectic  flight-tested  EKF  used  on  current  and  pre¬ 
vious  SSTL  missions14' 15  were  made. 

SATELLITE  ATTITUDE  PROBLEM 

Fundamentally,  the  problem  of  three-axis  attitude  parameterization  is  to  specify  the 
orientation  of  the  satellite  body  frame  with  respect  to  another  known  coordinate 
frame  (commonly,  the  local  orbital  or  inertial  coordinate  frames) .  Euler  angles  offer 
an  intuitive  way  to  represent  a  satellite’s  attitude.  These  angles  are  determined 
from  a  series  of  positive  right-hand  rotations  from  the  local  orbital  coordinate  frame 
to  the  satellite  body  frame.  In  this  paper,  the  pitch-roll-yaw  (2-1-3)  sequence  is 
used.  However,  though  Euler  angles  are  particularly  useful  due  to  their  conve¬ 
nient  physical  interpretation,  there  are  no  all-attitude  3-parameter  sets  that  are  free 
of  singularities.16  As  a  result,  it  is  often  convenient  to  utilize  the  singularity-free 
Euler  parameters  in  attitude  equations,  represented  in  this  paper  by  the  4-parameter 
quaternion  q. 

DYNAMICS 

Satellite  dynamics  in  inertial  space  are  governed  by  Euler’s  equations  of  motion. 
Describing  the  rotational  motion  of  the  satellite  with  respect  to  an  inertial  frame 
gives17 

IoJb/i  =  -utyi  x  (Iub/i  +  h)  -  h  +  N  (13) 

where 


In 

h  2 

hz 

hi 

I22 

hz 

hi 

I'12 

hz 

is  the  moment  of  inertia  tensor  of  the  satellite, 

ktyi  =  [wi  U>3]T  (15) 

is  the  inertially  referenced  body  angular  rate  vector, 

h  =  [hr  h2  h3]T  (16) 

is  the  satellite’s  internal  angular  momentum  vector  (produced  by  reaction  wheels 
or  control-moment  gyroscopes,  for  example),  N  is  the  total  external  torque  vector 
about  the  satellite  center  of  mass  with  respect  to  the  inertial  coordinate  frame,  and 
subscripts  1,  2,  and  3  are  the  axes  of  the  satellite  body  coordinate  frame.  Assuming 
the  principal  moments  of  inertia  axes  are  along  the  body  axes,  the  off-diagonal  terms 
of  Equation  (14)  are  zero. 
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KINEMATICS 


Kinematics,  or  the  study  of  motion  independent  of  forces,  are  governed  by  a  sep¬ 
arate  equation  based  on  quaternions.  Assuming  a  satellite  in  a  circular  orbit,  the 
quaternion  evolution  is18 

q  =  i[Q(0l)-fr(o2)]q  (17) 


where  Oi  = 


P 


,  o2  =  [0  n  0 


0  q  -  (p)  q  = 


=  n*  (q)  P  = 


0]T,  n  is  the  orbital  mean  motion,  and 


Pa  Pi  -Pi  Pi 

’  9i  " 

~Pi  Pi  Pi  P2 

92 

P2  ~Pi  Pi  Pi 

93 

~Pl  ~P2  -Pi  Pi  _ 

.  . 

<?4  -<?3  92  Ql 

’  Pi  " 

Qi  94  -Ql  92 

P2 

-92  Qi  94  93 

Pi 

co 

1 

CN 

Cr< 

1 

s< 

1 

Pi 

(18) 


SYMPLECTIC  PROPAGATOR  FOR  SATELLITE  ATTITUDE 

Consider  the  case  of  a  triaxial  satellite  in  a  circular  orbit  with  mean  motion  n.  Its 
rotational  dynamics  are  governed  by  Equation  (13)  and  its  kinematics  are  governed 
by  Equation  (17).  For  the  sake  of  convenience,  to  is  taken  to  mean  u >b/i  (inertially 
referenced  body  angular  rates)  and  motion  is  described  with  respect  to  the  principal 
satellite  axes  Ji,  /2,  and  /3.  Assuming  no  external  or  internal  torques  (i.e.  N  =  0 
and  h  =  0),  then  the  system  is  Hamiltonian  and  two  conserved  values  of  motion  are 
the  rotational  energy 

1  3 

E=  (19) 

^  fc=i 

and  the  angular  momentum  of  the  satellite  in  the  inertial  frame 

L  =  D{nt)AT{Iu>  +  h)  (20) 

where  A  is  the  direction  cosine  matrix  (DCM)  from  the  local  orbital  frame  to  the 
body  frame,  and  D  is  the  DCM  from  the  inertial  frame  to  the  local  orbital  frame. 
The  symplectic  IMPR  method  will  exactly  conserve  these  quadratic  first  integrals 
over  time.2  Applying  the  second-order  IMPR  method  from  Equation  (2)  to  the 
satellite  dynamics  and  using  superscripts  to  denote  propagation  steps  gives 

^  +  A t  (aicJ2cJ3  +  H12UJ2  +  tf13cJ3  +  u{)  (21a) 

^  =  m 2  ^  T  A t  (a2cJiaJ3  +  i72 1^1  +  H23103  +  u2)  (21b) 

=  c 0^  +  A t  (a3cJicJ2  +  H31U1  +  H321A2  T  ^3)  (21c) 
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where 


h 


uk  =  ^yL,  k  =  1,2,3 
h 


a  1 


uk 

h- 


wlj)  +  OJ 


O+i) 


k  =  1,2,3 


h- 


■>  <+2 

A 


/3-/l 


"  1  ^3 

hi 


/1-/2 


ho 


TT  _  O  TT  _  rr  _  '“O  r_T  _  1  TT  _  *  TT  _ 

#  12  —  7-,  #13  —  7-,  #21  —  7-,  #23  — - 7-»  #31  — - 7~,  #32  — 

■<1  11  12  12  13 


hi 


(22) 

(23) 

(24) 

(25) 


At  is  the  timestep,  and  superscripts  (j)  and  (j  + 1)  indicate  propagation  steps.  Once 
the  dynamics  are  propagated  to  give  the  satellite  angular  velocity,  the  next  step  is 
to  propagate  the  kinematics  to  determine  the  satellite  attitude.  Assuming  that  the 
angular  velocity  vector  in  body  coordinates  is  constant  over  a  single  timestep,  a 
closed-form  solution  to  Equation  (17)  exists19 


q,-+1  =  (26) 

This  solution  can  be  treated  via  a  two  step  method6 

qs  =  e5Atn*(°2)q.  =  cos  qj  +  Sm^2  ^  Q*  (o2)  (27a) 

•  f  IcjIAt'N 

qj'+i  =  e5Aifl(0l)qs  =  cos  q*  + - i^i  (0l)  qs  (27b) 

The  question  now  arises  as  to  the  relationship  between  the  dynamics  and  the  kine¬ 
matics.  In  the  case  of  a  Hamiltonian  system  which  can  be  separated  naturally 

H(q,p)  =  HA(q,p)  +  HB  (q,  p)  (28) 


then  different  rules  can  be  applied  to  each  variable  subset.  This  is  the  approach 
used  in  partitioned  Runge-Kutta  methods  such  as  the  second-order  Stormer-Verlet 
(leapfrog)  method2,20 


,  (  At  dHA 
2  ’  dp 

dHB 


p=Pjj 


Pj+ 1  =  Si  Pj ,  At, 


dq 


_  ,  (  A  t  dHA 
Qj+ 1  ^  1^+5’  2  ’  dp 


p=pj+ 1/ 


(29a) 

(29b) 

(29c) 


where  fi  and  /2  refer  to  the  specific  method  used  to  propagate  each  variable.  Com¬ 
bining  this  leapfrog  method  with  the  propagation  equations  outlined  above  for  the 
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dynamics  and  kinematics  gives  a  composite  numerical  method  for  the  satellite  at¬ 
titude,  and  every  composite  method  is  symplectic  when  applied  to  Hamiltonian 
equations.21 

A  number  of  other  modifications  to  this  model  have  been  made  but  have  not  yet 
been  included  in  this  analysis.  Though  the  IMPR  discretization  for  the  dynamics 
preserves  the  true  phase  space  trajectory,  a  shift  in  phase  along  the  trajectory  is 
possible  and  a  time  transformation  has  been  derived  that  accounts  for  the  shifts.13 
Additionally,  gravity  gradient  torques  have  been  incorporated  into  the  dynamical 
equations.  Though  angular  momentum  and  rotational  energy  are  no  longer  constants 
of  motion,  a  Jacobi  constant  arises  and  it  is  preserved  by  the  symplectic  attitude 
propagator.  In  the  future,  other  torques  may  be  included  in  the  dynamical  model 
as  perturbations,  which  has  been  done  previously  in  symplectic  orbit  propagators.21 

6-STATE  KALMAN  FILTER 

The  symplectic  propagator  described  above  was  added  to  the  flight-tested  SSTL 
6-stat.e  Kalman  filter  algorithm  to  determine  the  relative  merits  of  the  SKF  over  the 
EKF  when  applied  to  the  satellite  attitude  problem.  However,  the  use  of  quater¬ 
nions  for  kinematic  parameterization  in  the  symplectic  propagator  presents  a  unique 
problem.  Though  quaternions  are  4-parameter  sets,  they  are  not  independent  due 
to  the  unit  norm  constraint.  It  is  therefore  possible  to  estimate  all  four  quaternion 
values,  or  to  estimate  three  quaternion  values  and  derive  the  fourth.  Both  meth¬ 
ods  have  merits,  but  the  three  quaternion  approach,  known  as  the  6-state  filter,  is 
generally  more  computationally  efficient.22,23  The  algorithm  outlined  below  closely 
follows  the  SSTL  derivation  by  Hashida.24 

Assuming  that  the  true  state  is  a  small  deviation  from  the  predicted  state,  error 
values  can  be  defined  based  on  these  relationships 

q  =  5q  <g>  q_  (30a) 

up  =  5u>  +  uj~  (30b) 

where  the  *  subscript  indicates  a  state  estimate  and  the  -  subscript  indicates  that 
the  estimate  has  been  propagated  forward  in  time  to  give  a  state  prediction.  In 
addition  to  the  true  7-parameter  state,  x  =  [qr  uT]7  ,  the  6-parameter  error  state 
is  defined  as 

x  =  [5^  Ju;r]T  (31) 

where  the  ~  superscript  denotes  variables  associated  with  the  error  state;  5q  is  de¬ 
fined  by  5q  —  [5qT  5f/4]T;  and  Sq4  =  +yj  1  —  \Sq\2.  Linearizing  the  dynamics  and 
kinematics  from  Equations  (13)  and  (17)  and  using  the  same  assumptions  from  the 
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symplectic  propagator  gives 


5^=  -C{Cj)5q+huj  (32a) 

ISu  =  -[C(w)I-C  ( ICj  +  h)]  5u  (32b) 


where 


axb  =  C(a)b,  C{ a)  = 


0  — a3  a2 
o3  0  —  cq 
— a2  ai  0 


(33) 


These  linearized  equations  are  then  used  to  form  the  state  transition  matrix  for  the 
error  state 


C  (^n)  2^3x3 

03x3  -I-l[C{Cjn)I-C{I  u>n  +  hn)] 

(34) 

In  this  implementation,  direct  observation  of  <5<f  is  assumed  for  simplification,  giving 
the  following  observation  matrix 


—  16x6  T-Tnz\t,  Pn  — 


d 


^9, J 


d  (5q,  5u>) 


H  —  [I3x3  03x3] 


(35) 


The  6-state  Kalman  Filter  algorithm  follows  Equations  (4)  to  (8),  with  some  modi¬ 
fications: 

•  Calculate  the  current  time’s  7-parameter  state  prediction  x~  =  [q“T  u>'r]T 


rtn 

X“=Xn-l+  /  /  (xn-l)  dt 

Jtn- 1 


(36) 


•  Propagate  the  error  state’s  error  covariance  matrix  forward  in  time  using  the 
error  state  transition  matrix  in  Equation  (34) 

P;=*nPn-l*Z  +  Qn  (37) 

•  Calculate  the  error  Kalman  gain 

kn  =  P~Hl  [Hnp-Hl  +  _1  (38) 


•  Calculate  the  error  state  estimate  x„ 


x„  =  Kn  (zn  -  h  (x„)) 


(39) 


•  Calculate  the  error  state’s  error  covariance  matrix  estimate 

Pn  =  (l6x6  -  KnHn)  P-  (40) 
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•  Calculate  the  7-parameter  state  estimate  xn  =  [qj,  a>J] T  using  the  error  state 
estimate  xn,  the  7- parameter  state  prediction  x“,  and  Equations  (30) 

(41a) 

(41b) 

(41c) 
(41d) 

The  standard  6-state  EKF  uses  a  second-order  nonsymplectic  Adams-Bashford  method 
to  propagate  the  state  in  Equation  (36),  while  the  SKF  uses  the  symplectic  propa¬ 
gator  outlined  previously. 

SYMPLECTIC  KALMAN  FILTER  RESULTS 

To  compare  the  two  estimation  methods,  realistic  satellite  configurations  and  mis¬ 
sion  profiles  were  selected  for  simulation.  Many  configurations  were  tested,  and  a 
representative  case  was  based  on  the  upcoming  U.S.  Air  Force  Academy  FalconSAT- 
3  small  satellite  mission,25  with  one  modification.  In  order  to  test  the  findings  of 
Section  ,  a  high-accuracy  measurement  device  was  substituted  for  the  standard  low- 
accuracy  magnetometer  planned  for  FalconSAT-3.  Instead,  the  Advanced  Stellar 
Compass  (ASC)26  recently  flown  on  ESA’s  PROBA-1  mission  was  selected  due  to 
its  successful  implementation  as  a  rate-gyro  replacement  onboard  an  agile  resource- 
constrained  small  satellite,  and  for  its  ability  to  provide  high-accuracy  attitude 
observations  with  high  frequency  at  fast  rotation  rates  (up  to  2.5  (~^).  Table  3 
summarizes  the  simulated  satellite  properties. 


I  (kgm2) 

Mass  (kg) 

Obs  Freq  (Hz) 

Obs  Accuracy  (1  a,  arcsec) 

3.89  0  0 

0  3.89  0 

0  0  1.32 

50 

2 

1 

Table  3:  SKF  Versus  EKF  Simulation  Satellite  Properties 


In  addition,  properties  for  a  circular  orbit  were  derived  from  the  FalconSAT-3 
mission27,28  and  are  summarized  in  Table  4. 


Mean  Motion  (^) 

Orbit  Period  (min) 

Altitude  (km) 

Inclination  (deg) 

0.0627 

95.7 

560 

35 

Table  4:  SKF  Versus  EKF  Simulation  Orbit  Properties 


An  SSTL  RK-45  integrator  was  used  to  propagate  the  satellite  kinematic  and 
dynamic  equations  to  obtain  both  the  true  trajectory  and  the  observation  set  using 
observation  properties  from  Table  3  and  initial  conditions  from  Table  5. 


4 -+yi- 
<5q «  =  fe  Sq. 


Sqn  I 

lT 


qn  =  ®  q„ 

Cjn  =  5un  -F  Cj~ 
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9i 

92 

93 

94 

W\ 

B1 

w2 

(1H) 

\  sec  ^ 

: 

ws 

('deg 'i 

V  sec 

: 

At  (sec) 

0.0088 

0.0086 

0.0086 

0.9999 

2 

0.3 

0.5 

0.01 

Table  5:  SKF  Versus  EKF  Simulation  Initial  Conditions 

For  filter  initialization,  the  diagonal  values  of  the  covariance  matrix  for  both  the 
EKF  and  SKF  were  set  high  relative  to  the  measurement  noise  since  only  steady  state 
performance  was  considered.  Measurement  noise  was  selected  to  give  optimal  filter 
performance  and  process  noise  was  set  to  zero  to  minimize  masking  of  linearization 
errors,  consistent  with  previous  authors’  analysis.8, 17  The  state  vector  for  both 
filters  was  initialized  to  zero  values,  and  a  standard  small  satellite  ADCS  step  size 


of  0.1  seconds  was  used.  In  Figure  3,  the  significant  difference  between  the  two 


filters’  state  estimation  performances  can 


Time  (sec) 


seen. 


0 

2000 

4000  6000  6000 

Time  (sec) 

10000 

.  x  10"* 

EKF  Angular  Rate  Error 

Time  (sec) 


Figure  3:  SKF  Versus  EKF  State  Estimation  Error 


The  contrast  in  kinematic  and  dynamic  estimation  performance  is  highlighted 
using  a  logarithmic  scale  for  the  error  values  in  Figure  4.  Remarkably,  the  SKF 
state  error  is  one  order  of  magnitude  smaller  than  the  EKF  state  error  in  the  steady 
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state. 


Log  of  the  Absolute  EKF  and  SKF  Quaternion  Errors 


-  EKF  Quaternion  Errors 

—  SKF  Qualemion  Errors 


Time  (sec) 


Figure  4:  SKF  Versus  EKF  Log  State  Estimation  Error 


As  hypothesized,  the  SKF  preserves  constants  of  motion  to  a  much  higher  degree 
of  accuracy  than  the  EKF.  The  SKF  produced  rotational  energy  error  two  orders 
of  magnitude  smaller  than  the  EKF,  and  angular  momentum  errors  one  order  of 
magnitude  smaller  than  the  EKF.  These  results  are  presented  in  Figure  5. 
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Time  (sec)  Time  (sec) 


Figure  5:  SKF  Versus  EKF  Conserved  Values  Error 

The  SKF’s  better  steady-state  accuracy  and  conservation  of  constants  of  motion 
were  observed  throughout  the  simulations  conducted.  Only  by  reducing  the  EKF 
timestep  by  an  order  of  magnitude  in  this  case  was  it  able  to  match  the  steady- 
state  performance  of  the  SKF.  In  general,  the  filters’  performances  converged  as  the 
number  of  timesteps  per  observation  was  increased. 

CONCLUSIONS 

A  novel  method  for  efficient  high-accuracy  satellite  attitude  estimation  based  on 
symplectic  numerical  methods  has  been  presented.  The  application  of  a  symplectic 
integration  method  to  the  EKF  algorithm  outperformed  the  standard  EKF  in  the 
presence  of  nonlinearity  and  low  measurement  noise  in  the  1-D  case.  Using  this 
result,  a  six-state  SKF  was  developed  for  satellite  attitude  estimation.  Numerous 
small  satellite  simulations  were  conducted  and  the  SKF  demonstrated  orders  of  mag¬ 
nitude  improvement  in  state  accuracy  and  preservation  of  constants  of  motion  over 
the  EKF.  Based  on  these  results,  there  appears  to  be  a  clear  advantage  to  using  the 
SKF  in  place  of  the  EKF  onboard  small  satellites  that  do  not  have  the  computa- 
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tional  capacity  to  accommodate  increased  model  evaluations,  especially  those  with 
properties  modeled  above  (namely,  gyroless  agile  small  satellites  with  high-accuracy 
attitude  observation  devices). 
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Abstract: 

The  purpose  of  the  United  States  Air  Force  Academy  (USAFA)  Space  Systems 
Research  Program  is  to  give  cadets  the  opportunity  to  “learn  space  by  doing  space”  while 
also  providing  an  orbiting  platform  for  Air  Force  and  Department  of  Defense  (DoD) 
science  experiments  as  FalconSAT-3  is  designed  to  do.  This  paper  describes  small 
satellite  programs  at  the  U.S.  Air  Force  Academy’s  Space  Systems  Research  Center. 

FalconSAT-  2  and  Fa!conSAT-3  are  student-built  small  satellites  that  provide  low-cost  access  to 
space  for  DoD  space  research  &  development  payloads  as  well  as  platforms  for 
student  experiments.  Rapid,  low-cost  design  is  achieved  by  leveraging  Commercial  Off-the-  Shelf 
(COTS)  hardware  to  the  greatest  extent  possible.  FalconsSAT-2,  still  searching  for  an  alternate 
launch  opportunity,  was  the  first  to  demonstrate  the  use  of  COTS  modules  for  this  use. 
FalconSAT-3,  scheduled  to  be  launched  in  2006,  recently  completed  critical  design  and  built  upon 
the  successful  FalconSAT-2  experience  to  develop  an  even  more  capable  spacecraft  bus.  By 
using  the  off-the-shelf  equipment,  student  involvement  in  satellite  and  mission  design  has  been 
accelerated  and  provided  the  capability  to  challenge  students  This  paper  is  declared  a  work  of  the 
U.S.  Government  and  is  not  subject  to  copyright  protection  in  the  United  States,  through  more 
intense  participation  over  the  years.  Realizing  the  rapid  turnover  and  extended  commitments  of 
students  in  a  senior  undergraduate  program,  there  is  a  delicate  balance  to  be  found;  one 
between  comprehensive  mission  and  satellite  design  requirements,  and  adequate  experience  in 
a  multi-million  dollar,  real  world  space  program.  The  development  of  the  FalconSAT  program  will 
first  be  described  in  the  context  of  the  progress  made,  followed  by  a  more  detailed  discussion  of 
the  COTS  hardware  for  more  efficient  development  of  small  student  satellites  as  simple  payload 
platforms  for  educational  and  technological  purposes. 
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